#CHAPTER 5
##########

# simulating a DNA sequence

pi <- c(0.25,0.75)
x <- c(1,0)
set.seed(2009)
sample(x,10,replace=TRUE,pi)

pi<-c(0.25,0.25,0.25,0.25)
x<-c(1,2,3,4)
seq<-sample(x,10000,replace=TRUE,pi)
seq[1:15]

########

# simulating binomial random variables

x<- rbinom(2000,1000,0.25)
mean(x)
sd(x)^2
hist(x,xlab="Number of successes",main="")

x<- rbinom(10000,1000,0.25)
sum(x>=300)/10000
pbinom(300,size=1000,prob=0.25,lower.tail=FALSE)

bin25<-rbinom(1000,25,0.25)
av.bin25 <- 25*0.25
stdev.bin25 <- sqrt(25*0.25*0.75)
bin25<-(bin25-av.bin25)/stdev.bin25
hist(bin25,xlim=c(-4,4),ylim=c(0.0,0.4),prob=TRUE,xlab="Sample size 25",main="")
x<-seq(-4,4,0.1)
lines(x,dnorm(x))

pnorm(3.651501,lower.tail=FALSE)

qchisq(0.95,1)

########

# simulating a string according to Markov1 model

set.seed(2009)

markov1 <- function(x,pi,P,n){
mg <- rep(0,n)
mg[1] <- sample(x,1,replace=TRUE,pi)
for (k in 1:(n-1)){
  mg[k+1] <- sample(x,1,replace=TRUE,P[mg[k],])
  }
return(mg)
}

x<-c(1:4)
pi <- c(0.342,0.158, 0.158, 0.342)
# A C G T one-step transition matrix:
P <- matrix(scan(),ncol=4,nrow=4,byrow=T)
0.423 0.151 0.168 0.258
0.399 0.184 0.063 0.354
0.314 0.189 0.176 0.321
0.258 0.138 0.187 0.415

tmp <- markov1(x,pi,P,50000)
A<- length(tmp[tmp[]==1])
C<- length(tmp[tmp[]==2])
G<- length(tmp[tmp[]==3])
T<- length(tmp[tmp[]==4])
(C+G)/(A+C+G+T)   # fraction of G+C

count <-0
for (i in 1:49999){
  if (tmp[i]==2 && tmp[i+1]==3)
  count <- count+1
  }
count/49999 # abundance of CG dinucleotide as estimated by the model

########

# Poisson approximation to binomial

ppois(2,2.4)

########

# simulating restriction fragment lengths

set.seed(2009)

x<-c(1:4)
propn <- c(0.25,0.25,0.25,0.25)
seq2 <- sample(x,48500,replace=TRUE,prob=propn)
seq2[1:15]
length(seq2[])

rsite <- function(inseq, seq){
  # inseq: vector containing input DNA sequence,
  # A=1, C=2, G=3, T=4
  # seq: vector for the restriction site, length m
  # Make/initialize vector to hold site positions found in inseq
  xxx <- rep(0,length(inseq))
  m <-length(seq)
  # To record whether position of inseq matches seq
  truth <- rep(0,m)
    # Check each position to see if a site starts there
  for (i in 1:(length(inseq) - (length(seq) -1))){
    for (j in 1:m){
      if (inseq[i+j-1]==seq[j]){
      truth[j] <- 1 # Record match to jth position
      }
    }
    if (sum(truth[]) ==m){ # Check whether all positions match
    xxx[i] <- i          # Record site if all positions match
    }
    truth <- rep(0,m)    # Reinitialize for next loop cycle
  }
  # Write vector of restriction sites positions stored in xxx
  L <- xxx[xxx>0]
  return(L)
}

alu1 <- c(1,2,3,4)
alu.map <- rsite(seq2,alu1)
length(alu.map)
alu.map[1:10]

flengthr <- function(rmap,N){
  # rmap is a vector of restriction sites for a linear molecule
  # N is the length of the molecule
  frags <- rep(0,length(rmap))
  # Vector for substraction results: elements initialized to 0
  rmap <-c(rmap,N)
  # Adds length of molecule for calculation of end piece
  for(i in 1:(length(rmap)-1)){
  frags[i] <- rmap[i+1]-rmap[i]
  }
  frags <- c(rmap[1],frags) # First term is left end piece
  return(frags)
}
alu.frag <- flengthr(alu.map,48500)
alu.frag[1:10]
max(alu.frag[])
min(alu.frag[])
length(alu.frag[])
sum(alu.frag[])

hist(alu.frag,breaks=seq(0,2500,50), freq = TRUE,xlab="Fragment Size")

plot(c(0,2500),c(3,1),xlab="Fragment Size",ylab="",type="n",axes=F)
axis(1,c(0,500,1000,1500,2000,2500))
for (i in 1:length(alu.frag)){
lines(c(alu.frag[i],alu.frag[i]),c(1,3))
}

########

# k-word occurrences (basic version to perform the word count)

ec.prom <- read.table("promseqex1234.txt",sep="",header=F)
ec.prom <- as.matrix(ec.prom)
ec.prom <- ec.prom[,-ncol(ec.prom)]
ncol(ec.prom)
w <- 4 # restricting attention to 4-words

prob.ec <- c(0.246,0.254,0.254,0.246) # base frequencies for the E coli sequence
expect4.ec <- array(rep(0,4^w),rep(4,w)) # 4 is the max value in each dim for w
                                         # there are w dimensions
for (i in 1:4){
  for (j in 1:4){
    for (k in 1:4) {
      for (m in 1:4) {
        expect4.ec[i,j,k,m] <- 8*97*prob.ec[i]*prob.ec[j]*prob.ec[k]*prob.ec[m]
                            # 8 is the number of sequences in this example
                            # L-w+1 = 100 - 4 + 1 = 97
       }
    }
  }
}


Ncount4 <- function(seq,w){
  # w is length of word
  tcount <- array(rep(0,4^w),rep(4,w))
  # array[4 times 4 times 4 times 4] to hold word counts, elements set to zero
  N <- length(seq[1,]) # length of each sequence
  M <- length(seq[,1]) # number of sequences

  ##
  # Count total number of word occurrences
  for (j in 1:M){ # loopiing over sequences
    jcount <- array(rep(0,4^w),rep(4,w))
    # array to hold word counts for sequence j
    for (k in 1:(N-w+1)){ # looping over positions
    jcount[seq[j,k],seq[j,k+1],seq[j,k+2],seq[j,k+3]] <-
     jcount[seq[j,k],seq[j,k+1],seq[j,k+2],seq[j,k+3]] +1
     # adds 1 if word at k, k+1, k+2, k+3 appears in sequence j
    }
    tcount <- tcount + jcount
    # add contribution of j to total
  }
  return(tcount)
}

  prom.count <- Ncount4(ec.prom,4)
  sum(prom.count)
  max(prom.count)
 (1:256)[prom.count[,,,]>=10]

kwordseq <- NULL
for (i in 1:4){
  for (j in 1:4){
    for (k in 1:4) {
      for (m in 1:4) {
        kwordseq <- c(kwordseq,paste(i,j,k,m,sep=""))
      }
    }
  }
}
kwordseq[(1:256)[prom.count[,,,]>=10]]
prom.count[prom.count[,,,]>=10]
expect4.ec[(1:256)[prom.count[,,,]>=10]]

###


# k-word occurrences

ec.prom <- read.table("promseqex1234.txt",sep="",header=F)
ec.prom <- as.matrix(ec.prom)
ec.prom <- ec.prom[,-ncol(ec.prom)]
ncol(ec.prom)

prob.ec <- c(0.246,0.254,0.254,0.246) # base frequencies for the E coli sequence
expect4.ec <- array(rep(0,4^w),rep(4,w)) # 4 is the max value in each dim for w
                                         # there are w dimensions
for (i in 1:4){
  for (j in 1:4){
    for (k in 1:4) {
      for (m in 1:4) {
        expect4.ec[i,j,k,m] <- 8*97*prob.ec[i]*prob.ec[j]*prob.ec[k]*prob.ec[m]
                            # 8 is the number of sequences in this example
                            # L-w+1 = 100 - 4 + 1 = 97
       }
    }
  }
}


Ncount4b <- function(seq,w){
  # w is length of word
  tcount <- array(rep(0,4^w),rep(4,w))
  # array[4 times 4 times 4 times 4] to hold word counts, elements set to zero
  ncount <- array(rep(0,4^w),rep(4,w))
  # array[4 times 4 times 4 times 4] holds number of sequences with one or more of each k-word
  N <- length(seq[1,]) # length of each sequence
  M <- length(seq[,1]) # number of sequences

  ##
  # Count total number of word occurrences
  for (j in 1:M){ # loopiing over sequences
    jcount <- array(rep(0,4^w),rep(4,w))
    # array to hold word counts for sequence j
    for (k in 1:(N-w+1)){ # looping over positions
    jcount[seq[j,k],seq[j,k+1],seq[j,k+2],seq[j,k+3]] <-
     jcount[seq[j,k],seq[j,k+1],seq[j,k+2],seq[j,k+3]] +1
     # adds 1 if word at k, k+1, k+2, k+3 appears in sequence j
    }
    tcount <- tcount + jcount
    # add contribution of j to total

    ##
    # plug-in: add 1 to ncount if word occurs W= once in j
    for (k in 1:4){
      for (l in 1:4){
        for (m in 1:4) {
          for (n in 1:4) {
            if (jcount[k,l,m,n]!=0){
           ncount[k,l,m,n] <- ncount[k,l,m,n]+1
           }
          }
        }
      }
    }

  ##
  }
  return(tcount,ncount)
}
  prom.count <- Ncount4b(ec.prom,4)
  prom.count$ncount[1,2,1,1]
  max(prom.count$ncount)
  
ec.sim <- matrix(nrow=5000,ncol=51)
for (i in 1:5000){
  ec.sim[i,] <- sample(x,51,replace=T,prob.ec)
}
sim.count <- Ncount4b(ec.sim,4)
sim.count$ncount[1,1,1,1]/5000
sim.count$ncount[1,1,4,4]/5000
sim.count$ncount[2,1,2,4]/5000
#0.1356   0.1474  0.1642
prom.count$ncount[1,1,1,1]
prom.count$ncount[1,1,4,4]
prom.count$ncount[2,1,2,4]

Nw <- c(prom.count$ncount[1,1,1,1],prom.count$ncount[1,1,4,4],prom.count$ncount[2,1,2,4])
pw <- c(sim.count$ncount[1,1,1,1]/5000,sim.count$ncount[1,1,4,4]/5000,sim.count$ncount[2,1,2,4]/5000)
options(digits=4)
Zw <- (Nw-8*pw)/sqrt(8*pw*(1-pw))
1-pnorm(Zw)
